#include "SLAU.h"

#define line res->pointerB
#define col res->column
#define value res->val

CRSMatrix * allocation(unsigned int dx, unsigned int dy, unsigned int dz)
{
    unsigned int ddx=dx-2, ddy=dy-2, ddz=dz-2;
    CRSMatrix * res=new CRSMatrix;
    res->numRow=ddx*ddy*ddz;
    res->numColumn=res->numRow;
    line=new unsigned int[res->numRow+1];
    unsigned int count=1;

    line[0]=0;
    line[count++]=line[count-1]+4;              //line begin
    for(unsigned int k=2; k<ddx; k++)
    {
        line[count++]=line[count-1]+5;
    }
    line[count++]=line[count-1]+4;

    for(unsigned int j=2; j<ddy; j++)
    {
        line[count++]=line[count-1]+5;
        for(unsigned int k=2; k<ddx; k++)
        {
            line[count++]=line[count-1]+6;
        }
        line[count++]=line[count-1]+5;
    }

    line[count++]=line[count-1]+4;
    for(unsigned int k=2; k<ddx; k++)
    {
        line[count++]=line[count-1]+5;
    }
    line[count++]=line[count-1]+4;

    for(unsigned int i=2; i<ddz; i++)
    {
        line[count++]=line[count-1]+5;
        for(unsigned int k=2; k<ddx; k++)
        {
            line[count++]=line[count-1]+6;
        }
        line[count++]=line[count-1]+5;

        for(unsigned int j=2; j<ddy; j++)
        {
            line[count++]=line[count-1]+6;
            for(unsigned int k=2; k<ddx; k++)
            {
                line[count++]=line[count-1]+7;
            }
            line[count++]=line[count-1]+6;
        }

        line[count++]=line[count-1]+5;
        for(unsigned int k=2; k<ddx; k++)
        {
            line[count++]=line[count-1]+6;
        }
        line[count++]=line[count-1]+5;
    }

        line[count++]=line[count-1]+4;
    for(unsigned int k=2; k<ddx; k++)
    {
        line[count++]=line[count-1]+5;
    }
    line[count++]=line[count-1]+4;
    
    for(unsigned int j=2; j<ddy; j++)
    {
        line[count++]=line[count-1]+5;
        for(unsigned int k=2; k<ddx; k++)
        {
            line[count++]=line[count-1]+6;
        }
        line[count++]=line[count-1]+5;
    }

    line[count++]=line[count-1]+4;
    for(unsigned int k=2; k<ddx; k++)
    {
        line[count++]=line[count-1]+5;
    }
    line[count]=line[count-1]+4;                  //line end

    col=new unsigned int[line[count]];
    value=new T[line[count]];

    count=0;

    for(unsigned int i=0; i<res->numRow; i++)       //column begin
    {
        if(i/(ddx*ddy))
        {col[count++]=i-ddx*ddy;}
        if((i/(ddx))%ddy)
        {col[count++]=i-ddx;}
        if(i%ddx)
        {col[count++]=i-1;}

        col[count++]=i;

        if(i%ddx!=(ddx-1))
        {col[count++]=i+1;}
        if((i/(ddx)%ddy)!=(ddy-1))
        {col[count++]=i+ddx;}
        if(i/(ddx*ddy)!=(ddz-1))
        {col[count++]=i+ddx*ddy;}
    }
    return res;
}

void fillMatrix(CRSMatrix* A, T k, T ax, T ay, T az, T dt, T dx, T dy, T dz, unsigned int ddx, unsigned int ddy, unsigned int ddz)
{
    unsigned int count=0;
    ddx-=2; ddy-=2; ddz-=2;
    for(unsigned int i=0; i<A->numRow; i++)
    {
        if(i/(ddx*ddy))
        {
            A->val[count++]=-az*dt/(dz*dz);
        }
        if((i/(ddx))%ddy)
        {
            A->val[count++]=-ay*dt/(dy*dy);
        }
        if(i%ddx)
        {
            A->val[count++]=-ax*dt/(dx*dx);
        }

        A->val[count++]=1+k*dt+2*(ax*dt/(dx*dx)+ay*dt/(dy*dy)+az*dt/(dz*dz));

        if(i%ddx!=(ddx-1))
        {
            A->val[count++]=-ax*dt/(dx*dx);
        }
        if((i/(ddx)%ddy)!=(ddy-1))
        {
            A->val[count++]=-ay*dt/(dy*dy);
        }
        if(i/(ddx*ddy)!=(ddz-1))
        {
            A->val[count++]=-az*dt/(dz*dz);
        }
    }
}